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ABSTRACT 

The  problem  of  determining  the  singularity  expansion  method 
(SEM)  parameters  from  transient  thin-wire  data  is  examined.  A 
computer  code  is  used  to  generate  response  data.  The  SEM  param¬ 
eters  are  computed  from  these  data.  For  noisy  data,  it  is  shown 
that  the  parameter  values  can  be  improved  by  signal  averaging. 
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The  problem  of  determining  the  singularity  expansion 
method  (SEM)  parameters  from  transient  thin-wire  data  is  exam¬ 
ined.  A  computer  code  is  used  to  generate  response  data.  The 
SEM  parameters  are  computed  from  these  data.  For  noisy  data, 
it  is  shown  that  the  parameter  values  can  be  improved  by  signal 
averaging. 
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CHAPTER  I 


INTRODUCTION 

The  singularity  expansion  method  (SEM)  LI]  has  proven  to 
be  a  useful  tool  for  the  analysis  of  electromagnetic  scattering 
problems.  The  SEM  gives  a  relation  between  the  incident  field 
and  the  induced  currents  and  charges  on  a  body.  In  this  rela¬ 
tion,  the  natural  response  is  expressed  as  a  sum  of  complex  ex¬ 
ponentials.  The  SEM  parameters  in  the  sum  (natural  frequencies, 
natural  modes,  and  coupling  coefficients)  together  with  the  in¬ 
cident  waveform  describe  the  scattering  problem  in  a  useful  and 
compact  way. 

A  number  of  investigators  have  analyzed  simple  bodies  using 
SEM.  In  the  original  SEM  paper,  Baum  [2]  treated  a  perfectly  con¬ 
ducting  sphere  in  free  space.  He  obtained  analytical  expressions 
for  surface  current  density  and  surface  charge  density.  The  SEM 
parameters  were  found  numerically  from  these  expressions.  Tesche 
[3]  applied  the  method  of  moments  to  Pockl ington ' s  equation  and 
computed  the  SEM  parameters  for  a  thin  wire  in  free  space.  Wilton 
and  Umashankar  [4]  analyzed  an  L-shaped  wire.  Marin  [5]  derived 
approximate  expressions  for  the  natural  frequencies  and  natural 
modes  of  several  thin-wire  structures.  Marin  and  Liu  [6 J  used 
some  of  these  results  to  suggest  a  simple  way  of  solving  tran¬ 
sient  thin-wire  problems. 
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All  the  work  mentioned  so  far  is  theoretical  in  the  sense 
that  the  SEM  parameters  were  found  by  solving  equations  relating 
the  incident  field  and  the  corresponding  body  response.  For  peo¬ 
ple  involved  in  electromagnetic  pulse  ( EMP )  testing,  it  is  of  in¬ 
terest  to  know  whether  the  SEM  parameters  can  be  computed  from 
measured  test  data. 

Van  Blaricum  and  Mittra  [7]  began  work  in  this  area.  They 
used  Pronv's  method  [SJ  to  compute  the  natural  frequencies  of 
computer-generated  transient  data.  Their  results  gave  impetus  to 
a  number  of  other  studies.  Lager,  Hudson,  and  Poggio  [9]  applied 
Prony's  method  to  several  test  cases  and  provided  some  guidelines 
for  applying  the  method  [10].  There  are  some  problems  involved 
in  using  Prony's  method  on  noisy  data  [11].  The  performance  of 
the  method  can  be  improved  by  using  an  iterative  technique.  Such 
a  scheme  has  been  applied  successfully  to  actual  EMP  test  data 
[12].  At  this  point  the  problem  of  how  to  compute  the  natural  fre¬ 
quencies  of  a  body  from  transient  data  appears  to  be  fairly  well 
understood . 

The  next  line  of  work  involves  developing  techniques  for 
computing  the  other  two  parameters  (natural  modes  and  coupling 
coefficients)  from  real  data.  Pearson  and  Robertson  [13]  showed 
that  the  natural  modes  for  a  thin  wire  can  be  computed  from  tran¬ 
sient  data.  They  generated  current  responses  using  a  thin-wire 
computer  code.  The  exciting  function  was  a  Gaussian  voltage 
pulse  applied  across  a  gap  in  the  wire.  They  were  able  to  compute 


natural  frequencies  and  natural  modes  that  agree  well  with  those 
found  by  Tesche  [3]. 

The  work  reported  here  is  a  continuation  of  the  efforts 
mentioned  above.  The  aim  is  to  develop  techniques  for  estimating 
the  SEM  parameters  from  real  test  data.  A  general  approach  for 
estimating  the  parameters  consists  of  measuring  current  waveforms 
at  a  number  of  different  locations  on  a  test  body  for  various  an¬ 
gles  of  incidence  and  polarization.  The  parameters  are  estimated 
from  these  current  waveforms.  At  Dresent,  there  are  not  enough 
appropriate  test  data  to  use  in  the  calculations.  As  an  alterna¬ 
tive,  the  same  th in-wire  computer  code  adopted  in  several  of  the 
references  was  used  to  generate  transient  current  data.  The  thin 
wire  is  a  simple  scatterer  whose  SEM  parameters  are  well  known. 

It  provides  a  convenient  examole  for  testing  the  general  methods 
outlined  in  tnis  report. 

Our  study  produced  two  major  results:  The  first  is  that 
the  SEM  parameters,  including  the  natural  frequencies,  natural 
modes,  coupling  coefficients ,  and  normal  l'zation  factors  can  be 
calculated  numerically  from  transient  data  generated  with  broad 
and  narrow  excitations.  The  second  result  is  that  by  using  sig¬ 
nal  averaging  techniques,  the  SEM  parameters  can  be  calculated 
from  simulated,  noisy  EMP  test  data. 
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CHAPTER  II 


FORMULATION 

2.1  SEM  FOR  A  THIN-WIRE  SCATTERER 

Consider  a  thin-wire  scatterer  with  radius  a  and  length  L, 
assumed  to  lie  along  the  z  axis.  The  wire  is  struck  by  a  plane 
electromagnetic  wave.  The  direction  of  propagation  k  of  the  in¬ 
cident  field  makes  an  angle  a  with  the  normal  to  the  center  of  the 
scatterer  as  shown  in  Figure  1.  The  time  origin  is  taken  at  the 
point  where  the  incident  wavefront  peak  is  Dq  meters  from  the  mid¬ 
point  of  the  scatterer.  The  incident  field  E  is  assumed  to  lie 
in  the  plane  of  the  figure.  The  induced  current  on  the  scatterer 
for  any  a  and  0  £  z  <_  L  is  represented  in  [1]  as 

I(t,‘3,z)  =  £Ra(6,z)  exp(sat)u(t-tQ)  +  g(t,o,z)  (1) 

a 

where  t  is  the  turn-on  time,  the  s^  are  the  natural  frequencies 
or  poles  of  the  scatterer,  the  R  are  the  corresponding  residues 

Ot 

at  s^,  and  g(t,e,z)  is  the  forced  response.  The  summation  over 
■i  represents  the  natural  response  of  the  scatterer.  According 
to  SEM  the  residues  can  be  factored  as 

R  (9,2)  =  n-Xn>)i  (z)f(s  )  (2) 

a  a  a  a  x 
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■mi 


- J 


*■*  , 


t 


where 


f(sj 


-  the  normalization  factor 


-  normalized  coupling  coefficient 


-  normalized  natural  mode 


-  Laplace  transform  of  the  exciting  function 


As  a  result,  the  natural  response  can  be  expressed  in  terms 
of  SEM  parameters  as 

lR(1(!‘.z)  exp(s  ^t)u(t-to)  *  ^z)f(sj  exp(s  ^t)u(t-tQ) 

a 


Since  the  values  of  the  induced  current  are  always  real,  the 
poles  s  ^  shown  in  Equation  1  occur  in  conjugate  pairs  or  lie  on  the 
real  axis  of  the  complex  s-plane.  The  poles  have  negative  real 
parts  and  are  generally  assumed  to  be  simple  for  a  thin-wire  scat- 
terer.  The  pole  values  are  uniquely  determined  by  the  wire  radius 
and  length. 

The  natural  modes  describe  the  spatial  amplitude  variation 
in  the  current.  They  are  normalized  so  that  at  the  maximum  magni¬ 
tude  points  they  are  real  and  equal  to  one.  The  modes  are  inde¬ 
pendent  of  the  angle  of  incidence  e  but  are  a  function  of  the 
spatial  coordinate  z  and  the  scatterer  radius  and  length.  In  [6], 
the  authors  developed  an  approximate  expression  for  the  natural 
modes  of  a  thin  wire  as 

i  ^(z)  =  sinUiiz/L) ,  0  ^  z  ^  L.  \  =  1 ,2,3 .  (4) 
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The  coupling  coefficients  describe  how  the  incident  field 
couples  to  a  scatterer.  They  depend  on  the  angle  of  incidence 
•3  and  are  independent  of  2.  In  [6],  the  coupling  coefficients 
were  approximated  as 


f ( s  )nmaxn^  (e)  =  /  sin(airz/L)E1nc(z,s  )dz 
a  a  a  1  2  a 

0 


E^nc(z,s^)  is  the  Laplace  transform  of  the  incident  tan¬ 
gential  electric  field  along  the  scatterer  and  is  given  in  [3]  as 


Elnc(2,s)  =  f(s)cos(e)  exp(-S2  sin(<?)/c) 


where  c  is  the  velocity  of  light.  Hence,  Equation  5  becomes 


n  n  (e)  =  /  sin(air2/L)cos(t>)  exp(-s  2  sin(e)/c)dz  (7) 
a  a  0  01 

For  convenience,  the  coupling  coefficients  are  normalized 

so  that  at  the  maximum  magnitude  points  they  are  real  and  equal 

to  one.  The  normalization  factor  nmax  determines  the  proper  mag- 

01 

nitude  and  phase  in  Equations  2,  5,  and  7. 


2.2  THE  PROPERTIES  OF  THE  EXCITING  WAVEFORMS 

In  general,  the  induced  currents  on  the  scatterer  can  be 
viewed  as  the  transient  response  to  an  incident  EMP.  Two  incident 
waveforms  with  different  properties  were  employed  to  generate  the 
data  used  in  this  report.  They  are  the  Gaussian  pulse  and  the 
double  exponential  pulse. 
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(1 )  The  Gaussian  Pulse 

The  Gaussian  pulse  can  be  used  to  obtain  a  wide  bandwidth 
exciting  function.  It  is  given  by 


f -j  ( t )  =  exp(-t‘7o1') 


where  a  is  a  constant.  The  Laplace  transform  of  f,  (t)  is  given  by 


f^  (s )  =  o/ir  exp(o  s  /4) 


The  Gaussian  pulse  has  a  broad  frequency  spectrum  which 
resembles  that  of  the  delta  function,  except  that  the  width  of 
the  latter  is  boundless,  whereas  that  of  the  Gaussian  pulse  de¬ 
pends  on  the  spread  parameter  o.  The  bandwidth  of  the  pulse  in¬ 
creases  with  decreasing  o. 

Since  the  Gaussian  pulse  cannot  be  represented  by  a  low 
order  sum  of  complex  exponentials,  care  must  be  taken  when  ana¬ 
lyzing  the  wire  response  with  Prony's  method  [10].  One  can  wait 
until  the  exciting  Gaussian  waveform  has  effectively  crossed  the 
scatterer.  Alternately,  one  can  choose  a  small  enough  that  the 
pulse  bandwidth  is  greater  than  the  response  bandwidth.  In  this 
case,  to  a  good  approximation  the  natural  response  of  Equation  3 
is  the  scaled  impulse  response  and  g(t,o,z)  is  zero. 


(2)  The  Double  Exponential  Pulse 

The  double  exponential  pulse  is  given  by 


j 


f2(t)  =  Y(e 


00) 


-at  -6t  \ 

-  e  ) 


=  0, 


t  >0 
otherwi se 


where  8  and  a  are  the  constants  which  determine  the  rising  and  de¬ 
caying  rates  of  the  pulse.  The  field  intensity  y  is  assumed  to  be 
unity.  The  Laplace  transform  of  Equation  10  is 


Since  f^(t)  is  a  sum  of  complex  exponentials,  Prony's  meth¬ 
od  can  be  used  for  extracting  response  poles  from  the  source  driven 
region  where  the  exciting  waveform  is  still  present.  The  poles  due 
to  the  exciting  pulse  simply  show  up  in  the  list  of  poles  describ¬ 
ing  the  response. 


2.3  METHODS  FOR  CALCULATING  THE  SEM  PARAMETERS 

Two  techniques  were  used  to  calculate  the  natural  frequen¬ 
cies:  Prony's  method  [7],  [9],  [13]  and  the  iterative  prefilter 
method  [14].  Prony's  method  was  used  to  analyze  the  "noiseless” 
response  data  generated  by  the  thin-wire  code.  The  method  fits 
a  sum  of  complex  exponentials  to  a  given  current  response.  When 
the  current  actually  is  a  sum  of  complex  exponentials,  the  poles 
and  residues  computed  are  just  those  of  the  given  response.  Our 
experience  has  been  that  this  method  works  very  well  with  com¬ 
puter  generated  signals  where  the  noise  is  small.  With  this  kind 
of  data  it  is  easy  to  recognize  the  natural  frequencies  of  the 
scatterer  by  simply  increasing  the  number  of  exponentials  fitted. 
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As  the  order  increases,  some  pole  values  become  approximately  con¬ 
stant  to  several  decimal  places  while  others  do  not.  Numerical  er 
rors  in  the  computed  response  data  cause  some  of  the  fitted  pole 
values  to  continue  changing.  The  poles  that  stabilize  as  the  fit¬ 
ted  order  increases  are  the  natural  frequencies  of  the  thin-wire. 

In  actual  EMP  test  situations,  errors  occur  in  measuring 
the  response  of  an  object.  Any  practical  scheme  for  calculating 
the  SEM  parameters  must  be  able  to  tolerate  some  error  in  the  data 
There  are  a  number  of  error  sources  [15].  One  important  error 
arises  when  an  analog  waveform  is  digitized.  For  simplicity,  we 
will  neglect  other  errors  and  model  the  digitizing  error  as  a  se¬ 
quence  of  uncorrelated,  zero  mean,  random  variables  added  to  sam¬ 
ples  of  the  current.  Let  the  measured  samples  be  given  by 

y(n)  =  I(nAt,e,z)  +  e(n)  n  =  0,1,2, .  (12) 

where  I(nAt,9,z)  is  the  wire  current  at  sample  nAt  and  e(n)  is  the 
corresponding  error.  Prony's  method  can  be  used  to  estimate  the 
poles  of  I(nAt,e,z)  given  y(n).  But  it  is  necessary  to  use  a  num¬ 
ber  of  extra  poles  to  account  for  the  noise.  We  have  had  better 
results  using  the  iterative  prefilter  method  [14].  Briefly,  this 
technique  initializes  with  Prony's  method  and  then  iteratively 
filters  the  data  and  improves  the  difference  equation  coefficient 
(and  hence  the  pole)  estimates.  Our  experience  is  that  the  pre¬ 
filter  technique  gives  better  results  than  Prony's  method  when  the 
data  are  noisy. 
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The  calculation  of  the  natural  modes,  coupling  coefficients, 

and  their  normal i rat  ion  factors  depends  on  the  residues  R  (■>,;) 

1 

of  Equation  2.  These  residues  were  calculated  in  the  standard  wav 
[7]  after  the  poles  were  calculated  by  either  method  above.  It  is 
worth  mentioning  that  for  noisy  data  case  the  residue  estimates 
are  unbiased  if  the  poles  are  known  exactly  [16].  So  the  poles 
should  be  estimated  as  well  as  possible  before  the  residues  are 
calculated . 

The  data  records  I(t,;',;)  for  all  o  and  :  have  the  same 
poles.  So  there  is  a  redundancy  in  the  data  that  can  be  used  to 
advantage.  Also,  since  R^(o,r)  tactors  as  a  product  of  the  nat¬ 
ural  mode  and  coupling  coefficient,  there  is  a  redundancy  in  these 
parameters  as  well,  Section  3.2  in  Chapter  111  contains  a  dis¬ 
cussion  of  how  these  two  types  of  redundancy  can  be  used  to  im¬ 
prove  SEM  parameter  estimates.  The  averaging  schemes  described 
later  depend  on  the  assumption  that  error  in  the  originai  data 
has  zero  mean.  In  a  real  test  situation,  systematic  errors  can 
often  be  noticed  and  corrected,  so  this  assumption  is  reasonable. 
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CHAPTER  III 


NUMERICAL  RESULTS 

3.1  EVALUATION  OF  SEM  PARAMETERS  FROM  NOISELESS  DATA 

3.1.1  The  Transient  Response 

Since  sufficient  experimental  data  were  not  available  we 
simulated  data  by  running  the  time-domain  thin-wire  computer  code 
WT-MBA/LLL1B  [17].  The  thin-wire  scatterer  was  modeled  using  six¬ 
ty  equally-spaced  segments.  Two  different  exciting  fields  were 
used:  the  Gaussian  pulse  and  the  double  exponential  pulse.  The 
geometry  of  the  model  is  depicted  in  Figure  1.  For  both  cases, 
plane  wave  excitation  was  assumed;  DQ  is  the  distance  at  t  *  0 
between  the  peak  of  the  exciting  pulse  and  the  midpoint  of  the 
scatterer,  and  at  is  the  width  of  each  time  step.  These  param¬ 
eters  appear  in  the  thin-wire  code.  The  numerical  values  for  the 
parameters  of  the  models  are: 

(1)  The  Gaussian  Pulse  Excitation 

f^(t)  =  exp(-t2/c2) 

L  =  1  meter 
a  =  2  *  10"10  sec. 
a/L  =  0.005 

At  =  5.556  x  10  ^  sec. 

Dq  =  0.6  meter 

These  parameter  values  were  used  in  an  example  in  [7]. 


(2)  The  Double  Exponential  Pulse  Excitation 

f2(t)  =  (e*3t  -  e*et)u(t) 

L  =  10  meters 
»  =  4  x  106 
3  =  2.2  *  108 
a/L  =  0.005 

At  =  5.556  x  10‘10  sec. 

DQ  =  6  meters 

The  transient  response  referred  to  was  the  induced  current 
at  the  center  of  each  of  the  sixty  segments  and  was  computed  at 
500  time-steps  of  width  At.  The  SEM  expression  for  the  current 
is  given  by 

Kt.e.z)  =  In'Jaxn1(e)i  (z)f(s  )  exp(s  t)u(t-t  )  +  g(t ,e,z)  (13) 

&  a  oi  <x  a  o 

a 

Some  of  these  currents  have  been  plotted  in  Figures  2  through 
5,  showing  the  9  dependence.  In  Figure  6  through  9  the  effect  of 
the  variation  of  z  is  demonstrated.  The  time  origin  for  these 
graphs  was  taken  when  the  exciting  pulse  was  DQ  meters  away  from 
the  wire  center. 

The  damped  oscillatory  behavior  of  the  response  can  be  ob¬ 
served  from  all  of  the  previously  mentioned  figures.  The  responses 
to  the  Gaussian  pulse  have  more  high  frequency  content  than  the 
responses  to  the  double  exponential  pulse.  This  is  expected  from 
the  properties  of  the  exciting  functions.  The  response  to  the 
Gaussian  pulse  in  each  case  is  roughly  the  impulse  response  scaled 
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Figure  7.  Response  to  (a)  Gaussian  pulse;  (b)  double 
exponential  pulse  at  8  1  0°,  z/L  =0.5 


19 


Fiqure  8.  Response  to  (a)  Gaussian  pulse;  (b)  double 
exponential  pulse  at  •'  -  0° .  :,L  -  0.75 
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by  a  constant.  This  is  because  the  Gaussian  pulse  has  a  Fourier 
transform  that  is  more  or  less  flat  up  to  the  frequency  of  the 
seventh  wire  pole. 


The  responses  at  both  ends  of  the  scatterer  are  very  in¬ 
significant.  Theoretically,  the  responses  are  identically  zero 
since  the  responses  must  satisfy  the  boundary  conditions  for  the 
E-field  integral  equation  of  the  thin-wire  problem  [6]. 

I(t,*,0)  =  I(t.e.L)  =  0  (14) 

With  broadside  incidence  (6  =  0°),  the  electric  field  is 
parallel  to  the  wire  and  only  symmetric  modes  are  excited.  In 
this  case,  the  response  starts  exactly  at  the  same  time  ( =DQ/c 
sec.)  for  all  sixty  segments. 

Also,  as  the  incident  field  angle  approaches  90  degrees, 
the  response  tends  to  be  weakened.  This  is  because  the  exciting 
wave  travels  in  the  direction  of  z,  parallel  to  the  scatterer. 
Thus,  no  resonance  will  be  excited  on  the  scatterer  as  expected. 

3.1.2  The  Natural  Frequencies  and  the  Corresponding  Residues 

When  the  thin-wire  code  was  run  on  the  CDC  7600  computer, 
thirteen  decimal  precision  was  requested.  Hence,  the  output  data 
records  can  be  considered  practically  noise  free.  The  Prony  pro¬ 
gram  was  used  for  extracting  poles  and  calculating  the  residues 
successively  from  the  output  data. 
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Table  1  shows  the  natural  frequencies  obtained  by  running 


the  Prony  program  on  the  response  data  with  the  Gaussian  pulse 
excitation. 

Table  1 

First  Ten  Layer  One  Poles 

N  SGB^  Tesche^  Prony' s  Method 


1 

-0.0828 

+ 

j  0 . 9251 

-0.082 

+ 

jO. 926 

-0.0819 

+ 

J0.919 

2 

-0.1212 

+ 

j 1 .9117 

-0.120 

+ 

jl  .897 

-0.1210 

+ 

j 1  -  896 

3 

-0.1491 

+ 

j 2. 8835 

-0.147 

+ 

j  2 . 874 

-0.1490 

+ 

j  2 . 879 

4 

-0.1713 

+ 

j 3. 8741 

-0.169 

+ 

j  3 . 854 

-0.1717 

+ 

j  3 . 866 

5 

-0.1909 

+ 

J4.8536 

-0.188 

+ 

j4.835 

-0.1916 

+ 

j  4 . 852 

6 

-0.2080 

+ 

j 5. 8453 

-0.205 

+ 

j  5 . 8 1 7 

-0.2093 

+ 

J5.839 

7 

-0.2240 

•t- 

j 6. 8286 

-0.220  +  J6.800 

-0.2256 

+ 

j  6 . 82 1 

8 

-0.2383 

+ 

j  7 . 82 1 2 

-0.234 

+ 

J7.783 

-0.2386 

+ 

J7.794 

9 

-0.2522 

+ 

j 8. 8068 

-0.247 

+ 

J8.767 

-0.2575 

+ 

j  8 . 775 

10 

-0.2648 

+ 

j  9 . 800 1 

-0.260 

+ 

J9.752 

-0.2951 

+ 

j  9 . 7  35 

1 [ 1 8]  by  Singaraju,  Giri,  and  Baum 
2[3]  by  Tesche 

3with  the  Gaussian  pulse  at  9  =  60°  and  z/l  =  0.5 


The  results  reported  by  [3]  and  [18]  were  also  tabulated.  When  the 
first  ten  poles  were  compared,  our  results  agreed  very  well  with 
the  others.  These  pole  values  are  essentially  independent  of  the 
starting  time  for  the  analysis  and  also  of  the  segment  number  and 


angle  of  incidence  provided  the  corresponding  modes  are  excited. 
The  same  kind  of  results  are  obtained  when  the  double  exponential 
pulse  is  used  except  of  course  that  two  extra  poles  are  present. 

It  is  interesting  that  Tesche's  analysis  [3]  found  two  layers  of 
poles  but  our  approach  found  poles  from  only  the  first  layer. 

The  poles  and  residues  found  can  be  used  to  approximate  the 
original  waveform.  Figure  10  shows  the  reconstructed  responses  to 
both  excitations  at  0  degrees  incident  angle.  Nine  wire  pole 
pairs  were  used.  The  corresponding  transforms  are  shown  sepa¬ 
rately  in  Figures  11  and  12.  In  the  case  of  this  response  to  the 
double  exponential  pulse,  only  two  wire  pole  pairs  are  actually 
needed  to  provide  an  accurate  reconstruction  of  the  original  re¬ 
sponse.  Of  course,  the  number  of  pole  pairs  needed  in  the  recon¬ 
struction  is  a  function  of  the  frequency  content  of  the  original 
response.  This  frequency  content  depends  on  the  values  of  the 
natural  modes,  coupling  coefficients,  and  incident  field  spectrum. 

3.1.3  The  Natural  Modes 

For  a  fixed  e  =  e' ,  Equation  2  can  be  rewritten  as 

K  (e\z)  =  n^Vle'Ji  (z)f(s  ),  0  <  z  <  l  (15) 

a  a  a  a  a  —  — 

The  coupling  coefficient  term  in  Equation  15  is  merely  a 
constant  now.  As  z  is  varied  from  0  to  L,  R  (e'.z)  is  simply 
proportional  to  the  natural  mode  current  i z ) .  This  proportion¬ 
ality  may  be  removed  by  normalizing  the  maximum  value  of  each 
mode  to  unity. 
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Figure  12.  Fourier  transform  magnitude  of  reconstructed  response 
to  double  exponential  pulse  at  e  =  0°  and  z/L  =  0.5 


27 


The  normalized  natural  modes  for  the  first  three  resonances 
are  plotted  in  Figures  13  through  18.  It  is  observed  that  the  nat¬ 
ural  modes  are  either  symmetric  or  anti -symmetric  about  the  center 
of  the  scatterer.  When  the  field  is  normally  incident  on  the  scat¬ 
tered  only  even  modes  are  excited  as  expected. 

These  results  are  virtually  identical  with  those  obtained 
by  Tesche  [3]  and  by  Pearson  and  Robertson  [13]. 

3.1.4  The  Coupling  Coefficients 

As  can  be  seen  from  Equation  2,  the  coupling  coefficients 
contain  the  variation  of  the  residues.  For  a  fixed  z  we  calcu¬ 
lated  the  residues  for  several  different  The  value  of  z  can 
be  chosen  from  the  normalized  natural  mode  plot  so  that  the  peak 
of  each  mode  corresponds  to  the  chosen  position  of  z  on  the  scat¬ 
terer.  The  normalized  variation  in  these  residues  gives  the 

The  normalized  coupling  coefficients  for  the  first  three 
resonances  are  presented  in  Figures  19  and  21  for  both  excitations 
with  both  real  and  imaginary  parts  shown.  In  Figure  20,  the  imag¬ 
inary  part  (Gaussian  excitation)  is  shown  alone. 

It  should  be  noted  that  we  normalized  the  coupling  coef¬ 
ficients  so  that  the  maximum  magnitude  points  are  real  and  equal 
to  one.  This  is  convenient  when  handling  experimental  data  since 
the  underlying  integral  equation  is  unknown.  Tesche  normalized 
the  coupling  coefficients  differently  since  he  knew  the  integral 
equation.  Thus  different  versions  of  the  plots  were  obtained. 
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Imaginary  Darts  of  first  three  normalized  natural 
modes  for  thin  wire,  a/L  =  0.C05  (Gaussian  pulse 
excitation  at  e  =  60°) 
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Figure  15.  Real  parts  of  first  three  normalized  natural 
modes  for  thin  wire,  a/L  =  0.005  (double 
exponential  pulse  excitation  at  6  =  60°) 
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Figure  18.  Imaginary  parts  of  first  three  normalized  natural 
modes  of  thin  wire,  a/l  =  0.005  (Gaussian  pulse 
excitation  at  o  =  0°) 
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Figure  21.  Normalized  coupling  coefficients  calculated  using 
double  exponential  excitation,  both  real  and 
imaginary  parts  of  first  three  modes,  are  shown 
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As  can  be  observed  from  the  graphs,  the  imaginary  parts  of  the 
coupling  coefficients  are  relatively  small  and  this  implies  that 
the  coupling  coefficients  are  almost  real  functions  of  the  inci¬ 
dent  angle  0.  The  graphs  would  be  smoother  if  more  values  of 
n|(e)  had  been  computed.  Only  nineteen  values  of  0  were  used  for 
each  graph  in  Figures  19,  20  and  21. 


3.1.5  The  Normalization  Factors 


Now  all  the  SEM  parameters,  except  the  normalization  fac¬ 
tors,  in  Equation  2  are  known.  This  final  parameter  may  be  deter¬ 
mined  directly  from  Equation  2  by  means  of  simple  complex  arith¬ 
metic,  i  .e. , 


max 

1. 


T 

n 

a 


S  t, 

R  (e,z)e  J 
(e)i  (z)f(s  ) 

01  Ot 


for  0  <  z  <_  L  (16) 


The  normal ization  constant  contains  an  arbitrary  factor  t^ 
that  depends  on  the  choice  of  the  time  origin.  In  all  the  pre¬ 
ceding  graphs,  the  time  origin  was  chosen  so  that  at  t  =  0  the  in¬ 
cident  field  peak  was  DQ  meters  from  the  wire  center.  DQ  is 
arbitrary  and  is  chosen  for  convenience  in  running  the  thin-wire 
code.  However,  it  is  better  to  eliminate  D„  from  nmax  ond  have 
it  depend  on  the  shape  of  the  scatterer.  For  this  reason,  we 

multiplied  the  residues  R  (e,z)  by  a  time  shift  factor  to  move 

01 

the  time  origin  to  the  time  when  the  incident  pulse  peak  crosses 

Vl 

the  wire  center.  The  required  factor  is  e  ,  where 


ll  ■  Vc  -  *s 


(17) 


and  where  t  is  the  time  at  which  we  started  to  analyze  the  tran¬ 
sient  response.  In  our  calculations,  ts  was  chosen  at  the  210^ 
time  step  to  ensure  all  the  observed  responses  had  started.  So 
t  =  210  (At). 

The  normalization  factors  for  the  first  three  modes  of 
both  excitations  are: 

1.  The  Gaussian  Excitation 

n'jiax  =  -0.7072779*10^  -  jO.  19381 36*106  (0=0°  ,z/L=0. 5) 

n£3X  =  0. 1431957*106  -  j0.4662369*106  (0=60° ,z/L=0.25) 

n'J35*  =  0.3933103*106  +  j0. 1275519*106  (0=0° ,z/L=0.17) 

2.  The  Double  Exponential  Excitation 

nfJ1dX  =  -0.7075146*10^  -  jO.  1935304*106  (0=0°  ,z/L=0.5) 

n™3X  =  0.14211 89 x  1 0 6  -  j0.4659009*106  (0=60°  ,z/L=0.25) 

n™*  ■  0. 391 1 535*106  +  jO. 1245812*10^  (e»0°,z/L»0.17) 

Since  the  normalization  factor  is  merely  a  scaling  complex 
constant  for  each  pole,  it  is  independent  of  either  0  or  z.  The 
results  obtained  are  very  close  for  both  excitations. 

Tesche  [3]  normalized  the  coupling  coefficients  differently, 
therefore  his  normalization  factors  are  different  from  what  we 
computed. 
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3.2  THE  EFFECTS  OF  NOISE  ON  SEM  PARAMETER  DETERMINATION 


3.2.1  The  Effect  of  Noise  on  the  Transient  Responses 

The  simulated  EMP  data  were  obtained  by  adding  sequences 
from  the  CDC  6600  random  number  generator  to  the  data  computed  by 
the  thin-wire  code.  The  noise  sequence  was  assumed  to  be  uniformly 
distributed  with  (statistical)  mean  zero  and  standard  deviation 
scaled  to  one  and  one  half  percent  of  the  peak  value  of  each  data 
record.  This  is  about  the  level  of  digitizing  noise  encountered 
in  actual  testing. 

Sample  noise  corrupted  responses  (e=G°,z/L=0.5)  are  pre¬ 
sented  in  Figure  22  for  both  excitations.  The  corresponding  Fourier 
transform  magnitudes  are  shown  in  Figures  23  and  24. 

3oth  signals  contain  the  1.5  percent  standard  deviation 
noise.  More  signal  peaks  are  visible  in  the  Gaussian  response  spec¬ 
trum.  This  is  because  the  wide  bandwidth  Gaussian  pulse  puts  more 
energy  into  the  higher  frequency  modes  than  does  the  double  exponen¬ 
tial  pulse. 

3.2.2  The  Effect  of  Noise  on  the  Natural  Frequencies  3nd  the 

Residues 

The  direct  computation  of  the  poles  from  noisy  transient 
data  can  lead  to  inaccurate  results.  Besides,  only  3  limited  num¬ 
ber  of  system  poles  can  be  extracted,  since  the  rest  have  been 
severely  corrupted  by  noise.  In  principle,  each  response  record 
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Figure  24. 


Fourier  transform  magnitude  of  noise  corrupted 
response  to  double  exponential  pulse  at  0  =  0° 


z/L  =  0.5 


contains  the  same  poles.  It  is  possible  to  use  extra  data  to  en¬ 
hance  a  particular  natural  frequency. 

The  iterative  prefilter  method  rather  than  Prony's  method 
was  used  for  all  the  pole  estimates  in  this  section.  Suppose  data 
y(t,e,z)  is  available  at  several  e  and  z  values.  Poles  can  be  com¬ 
puted  from  each  of  the  records  resulting  in  estimates  s  ($,z).  The 
arguments  e  and  z  are  included  to  show  that  the  pole  estimates  de¬ 
pend  on  location  and  angle  of  incidence.  Of  course  the  actual  wire 
poles  are  independent  of  a  and  z.  Usually  some  of  the  pole  esti¬ 
mates  cluster,  that  is,  the  values  lie  close  together  in  the  s- 
plane.  An  estimate  of  s  can  be  computed  as  the  average  over  a 

CL 

cluster. 

=  ills  (e,z)  (18) 

a  n  ez  a 

where  N  is  the  total  number  of  poles  in  the  cluster.  We  initially 
planned  on  using  s  as  a  final  estimate  of  s  .  For  low  order,  say 
a  =  1,2,  the  results  are  adequate.  But  the  higher  order  poles  had 
enough  error  that  we  looked  for  a  way  to  improve  the  estimate. 

The  problem  with  using  s^  as  an  estimate  is  that  the  sa(e,z)  are 
biased  estimates  of  sa.  Averaging  tends  to  remove  the  random  com¬ 
ponent  of  error  but  can  do  nothing  about  the  bias.  It  would  be 
better  if  some  of  the  noise  could  be  averaged  out  of  the  original 
data  records  before  calculation  of  the  pole  estimates.  This  sug¬ 
gests  averaging  the  data  records  instead  of  the  s  (e,z).  The  data 

Cl 


records  used  in  the  averaging  ought  to  be  those  where  the  a -natural 
mode  is  near  a  peak.  That  is,  the  average  should  be  weighted  ap¬ 
proximately  in  proportion  to  the  natural  mode. 

As  a  next  step  we  used  the  as  poles  and  computed  residue 
estimates  R  (e,z)  from  the  data.  For  the  value  of  9  =  e  that  ap- 
peared  to  give  the  largest  range  of  variation  in  the  magnitude  of 
R  (e  ,z)  we  computed  an  estimate  of  the  normalized  natural  mode 

a 

i  (e„,z).  This  natural  mode  is  an  estimate  of  which  data  records 

a  0 

have  the  best  signal  to  noise  ratio  for  the  a-natural  frequency 
component.  The  data  y(t,6Q,z)  corresponding  to  z  values  where 
M60,z)  is  approximately  ±1  should  be  the  best  data  for  esti¬ 
mating  s  .  Let 

a 

*a(t)  3  l  Wa(z)y(t,eQ,z)  (19) 

be  a  weighted  average  of  the  data,  (j>a(t)  depends  on  0Q  but  we 
will  not  explicitly  show  this.)  The  sum  is  taken  over  the  avail¬ 
able  z  values.  One  choice  for  the  weight  function  W^(z)  is  to 
make  it  proportional  to  the  estimated  natural  mode  i  (eQ,z). 

This  proved  to  be  an  acceptable  choice  for  low  values  of  a  since 
the  corresponding  ia(90>z)  is  fairly  smooth.  But  with  a  =  5,  for 
example,  the  natural  mode  estimate  is  more  erratic.  As  a  result 
we  used  a  coarse  weight  function  that  was  proportional  to  +1  for 
z  values  that  make  i  (e„,z)  =  1,  equal  to  -1  for  z  values  that 
make  i aQ ,z )  =  -1  and  equal  to  zero  otherwise.  Figures  25,  26, 
and  27  show  graphs  of  the  $^(t)  obtained  for  a =  1,  3,  and  5.  The 
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weight  functions  Wa(z)  were  chosen  to  average  six,  nine,  and  fif¬ 
teen  data  records,  respectively.  The  graphs  indicate  that  the 
noise  level  increases  as  a  goes  from  1  to  5. 


The  pole  estimate  s  computed  from  c>  (t)  should  be  better 

at  a 

than  from  an  estimate  computed  from  any  one  of  the  y(t,a,z).  This 

is  because  the  noise  variance  is  reduced  by  the  average  over  z. 

Moreover,  the  pole  estimate  computed  from  $  (t)  should  be  better 

01 

than  an  estimate  computed  from  averaging  the  s  (o,z).  This  is  be¬ 
cause  the  s  (e,z)  contains  a  bias  error  while  the  original  data 

a 


records  y(t,6,z)  do  not  (according  to  our  assumption  about  the 
noise  having  zero  mean). 


As  an  example  of  the  results  from  using  the  -t> ^ ( t )  technique 
with  a  =  0°,  we  have 


s]  =  -0.0824553  +  jO. 9183957 
s3  =  -0.1581486  +  j2. 379185 
s5  =  -0.1893679  +  j4. 834777 


These  pole  estimates  still  have  some  error,  of  course,  but 
they  are  better  results  than  we  were  able  to  get  by  averaging  the 
s^(a,z)  values  computed  from  the  same  data. 

Figures  25,  26,  and  27  show  an  interesting  fact.  The  $  (t) 

a 

contain  mostly  the  a  natural  frequency.  To  see  why  this  is  the 
case,  consider  neglecting  g(t,e,z)  and  rewriting  Equation  3  as 

I(t,o,z)  =  ][*  (t)  s i n (cxttz/L ) ,  0  ^  z  ^  L 

a  a  =  1,2,3 .  (20) 
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where 


and  where  the  approximation  of  Equation  4  is  used  for  the  natural 
mode.  Equation  20  has  the  form  of  a  Fourier  series  in  z. 

The  Fourier  series  coefficients  ♦  (t)  are  given  by 

$  (t )  =  r  /  I(t,e,z)  sin(«7:z/L)dz  (22) 

*  1  o 

The  function  o  (t)  contains  time  behavior  corresponding  to 
01 

s  alone.  If  the  integral  in  Equation  22  is  approximated  by  a 
01 

sum  over  M  values  of  z  spaced  az  apart,  then  we  can  write 

O,.  M 

4  (t)  =  -J-  l  I(t,e,z)  sin(curZ./L)  (23) 

a  u  i=i  i 

where  the  zi  are  the  z  values.  This  equation  compares  with  Equa¬ 
tion  19.  The  functions  4  (t)  establish  a  bound  on  how  well  the 

01 

various  natural  frequencies  can  be  separated. 

The  4  (t)  were  computed  from  noisy  data  and  the  results  are 
01 

plotted  in  Figures  28,  29,  and  30  for  a  =  1 ,  3,  and  5.  In  these 
calculations  M  =  60  and  az  =  10/60  meter.  The  pole  values  s^ 
estimated  from  applying  the  iterative  prefilter  algorithm  to$Jt) 
are 

•  -0.08219866  +  jO. 9187515 
s  3  =  -0.1480808  +  j 2. 877230 
s5  =  -0.1933258  +  j4. 847712 
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As  would  be  expected,  these  results  are  better  than  those 
obtained  from  the  $  (t).  Unfortunately,  the  procedure  is  not 
practical  since  we  usually  do  not  know  the  natural  modes.  And  of 
course  the  modes  are  not  exactly  orthogonal.  The  thin  wire  is  a 
special  case  in  that  the  modes  are  approximately  orthogonal.  But 
the  principle  is  that  the  data  should  be  weighed  proportionally 
to  the  natural  modes. 

The  technique  for  calculation  of  s>  (t)  could  be  extended  by 
averaging  several  ^(t)  functions  computed  with  different  3  val¬ 
ues— perhaps  weighed  by  an  estimate  of  the  coupling  coefficient. 
This  was  not  tried. 

From  the  form  of  Equation  2,  it  is  clear  that  for  some 
values  of  8  and  z  not  all  the  modes  are  present  in  I(t,e,z). 

This  is  simply  because  of  zeros  in  n^(8)  and  i^(z) .  When  the 
data  are  noisy,  there  is  a  "dead  zone"  around  the  zeros  of  these 
functions  where  the  signal  to  noise  ratio  is  low,  and  the  natural 
frequencies  cannot  be  estimated  accurately.  For  experimental  data 
we  expect  that  it  may  be  necessary  to  use  a  two  step  measurement 
approach.  At  the  first  step  one  places  sensors  of  the  test  ob¬ 
ject  guided  by  some  pretest  calculations.  After  analyzing  the 
resulting  data  and  estimating  the  SEM  parameters,  the  second  step 
is  to  modify  the  sensor  placement  and  field  orientation  to  better 
identify  some  modes  or  to  find  missing  modes. 
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3.2.3  The  Effect  of  Noise  on  the  Natural  Modes 

In  order  to  prevent  the  occurrence  of  bias  error  resulting 
from  the  pole  estimation  process,  the  best  pole  values  obtained 
should  be  used  to  compute  the  residues  for  the  determination  of 
the  natural  modes.  If  the  pole  values  are  accurately  known,  then 
approximately  unbiased  solutions  for  residues  will  be  obtained. 

The  real  and  imaginary  parts  of  the  natural  modes  calcu¬ 
lated  by  using  the  pole  values  obtained  from  both  $  (t)  and 

Ol 

i  (t)  are  presented  in  Figures  31  through  34  for  a =  1,  3,  and  5. 

OL 

The  graphs  indicate  that  noise  is  more  of  a  problem  with  the 
higher  frequency  modes  than  with  the  fundamental. 

The  results  obtained  can  be  improved  by  averaging  the  nat¬ 
ural  modes  for  the  same  a,  computed  from  several  values  of  e. 

That  is,  one  can  compute  estimates  i  (e,z)  of  i  (z)  by  using  res- 

Cl  01 

idues  calculated  from  data  obtained  at  different  angles  of  inci¬ 
dence.  Then  an  estimate  of  the  natural  mode  is  obtained  as 

i  (z)  =  ff"  (e,z) 

N0  5  a 

where  N4  is  the  number  of  different  s  values  used.  This  estimate 
is  unbiased  if  the  residue  estimates  are  unbiased.  And  of  course 
averaging  reduces  the  error  variance.  The  average  could  include 
a  weight  factor  proportional  to  the  coupling  coefficient  but  this 
was  not  done. 
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The  averaged  natural  mode  for  a  =  3  is  shown  in  Figures  35 
and  36  for  real  and  imaginary  parts,  respectively.  Three  uncor¬ 


related  noisy  third  modes  were  averaged  to  produce  this  result. 

3.2.4  The  Effect  of  Noise  on  the  Coupling  Coefficients 

The  coupling  coefficients  were  computed  by  using  the  poles 
obtained  from  $  (t)  for  the  reason  given  in  the  previous  section. 
The  real  and  imaginary  parts  of  the  coupling  coefficients  for 
a  =  1,  2,  and  3  are  plotted  in  Figures  37  and  38,  respectively. 
The  coupling  coefficient  for  a  =  3  was  improved  by  averaging 
three  noisy  third  modes.  The  average  has  the  form 

■\(d)  =  f7-  [nVe.z)  (25) 

2  Z 

where  n^(o,z)  is  the  normalized  coupling  coefficient  at  location 
z,  and  N,  is  the  number  of  z  values  used. 

4. 

These  results  are  presented  in  Figure  39,  the  real  part; 
and  Figure  40,  the  imaginary  part. 

3.2.5  The  Effect  of  Noise  on  the  Normalization  Factors 

The  normalization  factors  can  still  be  determined  by  Equa¬ 
tion  16  since  both  the  natural  modes  and  coupling  coefficients 
are  now  known.  However,  experience  shows  that  because  of  resi¬ 
due  error,  a  more  accurate  value  will  be  achieved  by  averaging 
several  nmax  for  the  same  which  have  been  calculated  from  dif- 

Cl 

ferent  values  of  a  and  z.  Averaging  three  noisy  factors  for  each 


Figure  35.  Real  part  of  normalized  natural  modes  for 
.>  =  3,  obtained  by  signal  averaging 
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Figure  36.  Imaginary  part  of  normalized  natural  modes 
for  a  =  3,  obtained  by  signal  averaging 


case,  we  obtained  the  normal ization  factors  for  the  first  three 
resonant  frequencies  as 

n™ax  =  -0.71 48283*1 06  -  jO. 1935363*106 

n™*  =  0.1 383078*10®  -  j0.4762430*106 

n™X  =  0.4027178*10^  +  jO. 1 609751 *106 

Implicit  in  the  n™ax  calculation  is  the  assumption  that 
f(s  )  is  known  exactly.  This  assumption  is  unrealistic  for  the 
higher  frequency  poles.  In  EMP  testing,  information  on  f(t)  is 
obtained  from  field  map  data.  From  this  data  we  can  expect  to 
estimate  the  low  frequency  behavior  of  f(s)  fairly  well  but  not 
the  high  frequency  behavior.  So  how  accurately  f ( )  is  known 
depends  on  s^.  For  the  low  frequency  poles  of  a  large  object, 

we  expect  only  minor  problems  with  the  f(s  )  values. 

01 
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CHAPTER  IV 


CONCLUSIONS 

The  SEM  parameters  for  a  thin  wire  can  be  evaluated  di¬ 
rectly  from  noiseless  transient  data  using  Prony's  method.  Our 
results  agree  well  with  Tesche's  [3].  With  noisy  data,  the  iter¬ 
ative  prefilter  method  gives  better  results.  Averaging  tech¬ 
niques  can  be  used  to  enhance  the  values  of  the  parameters.  We 
tested  two  ways  to  improve  the  pole  values.  The  first  approach 
works  simply  by  averaging  several  pole  values  for  the  same  a 
which  have  been  computed  from  the  currents  at  different  loca¬ 
tions  on  the  wire  and  different  field  orientations.  This  algo¬ 
rithm  gave  reasonably  good  results  for  the  lower  frequency  poles. 
The  second  approach  involves  averaging  several  data  records  which 
have  been  weighed  approximately  in  proportion  to  the  natural 
modes.  The  poles  can  be  extracted  directly  from  these  averages. 
This  second  technique  gives  more  accurate  pole  values  than  the 
first. 

It  is  interesting  to  note  that  only  layer  one  poles  [3] 
were  found.  The  number  of  these  poles  needed  to  reconstruct  the 
natural  response  of  a  thin-wire  scatterer  depends  on  the  values 
of  the  natural  modes,  coupling  coefficients,  and  incident  field 
spectrum.  For  some  responses,  a  two  pole-pair  reconstruction 
is  adequate;  other  responses  require  ten  or  more  pole  pairs. 


The  best  pole  values  should  be  used  in  the  determination 
of  the  natural  modes  and  coupling  coefficients  to  minimize  the 
effect  of  bias  error.  Then  these  parameter  estimates  can  be 
improved  by  signal  averaging.  The  values  of  the  corresponding 
normalization  factors  can  be  improved  by  the  same  process.  When 
these  techniques  are  applied  to  an  arbitrary  test  object,  it  may 
be  necessary  to  use  a  two  (or  more)  step  process.  At  the  first 
step  sensor  placement  and  incident  field  orientation  are  chosen 
based  on  pretest  calculations  and  scale  model  results.  After 
estimating  the  SEM  parameters  from  the  resulting  data,  the  second 
step  is  to  modify  the  sensor  placement  and  field  orientation  to 
better  identify  some  modes  or  to  find  missing  modes. 

The  next  phase  of  our  research  will  be  to  treat  more  com¬ 
plicated  structures  which  provide  a  better  approximation  to  air¬ 
craft  configurations.  Our  ultimate  goal  will  be  to  see  how  the 
calculations  are  affected  when  coping  with  actual  aircraft  test 
data. 
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